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Abstract. The proper modeling of nonequilibrium gas dynamics is required in certain regimes of 
hypersonic flow. For inviscid flow this gives a system of conservation laws coupled with source terms 
representing the chemistry. Often a wide range of time scales is present in the problem, leading to numerical 
difficulties as in stiff systems of ordinary differential equations. Stability can be achieved by using implicit 
methods, but other numerical difficulties are observed. The behavior of typical numerical methods on a model 
advection equation with a parameter-dependent source term is studied. Two approaches to incorporate the 
source terms are utilized: MacCormack type predictor-corrector methods with flux limiters, and splitting 
methods in which the fluid dynamics and chemistry are handled in separate steps. Various comparisons over 
a wide range of parameter values are made. In the stiff case where the solution contains discontinuities, 
incorrect numerical propagation speeds are observed with all of the methods considered. This phenomenon 
is studied and explained. 


1. Introduction. In nonequilibrium gas dynamics, chemical reactions between the 
constituent gases must be modeled along with the fluid dynamics. This added complexity 
is required in certain regimes of hypersonic aerodynamic modeling, for example in the bow 
shock of hypersonic vehicles. 

Coupled systems of this form also arise in combustion problems. In particular, the mod- 
eling of scramjet engines that might be used in hypersonic vehicles requires the numerical 
simulation of supersonic combustion. 

Restricting our attention to inviscid flow, we have essentially the Euler equations of gas 
dynamics, coupled with source terms representing the chemistry. In two space dimensions 
these equations take the form 

(1) u t + /( u) x + g(u)y = ip(u) 

where u is the vector of dependent variables including momentum, energy, and densities or 
concentrations for each species in the reacting mixture. The flux functions / and g describe 
the fluid dynamics as in the Euler equations while the source term ip(u) arises from the 
chemistry of the reacting species. 

A variety of such systems are possible, depending on the level of detail of chemical 
modeling included. Examples and more discussions of these equations may be found in 
various references, e.g. [2], [13], [18]. 

When we attempt to solve the reacting flow equations numerically, new difficulties arise 
that are absent in non-reacting flows. Aside from the increase in the number of equations, 
the main difficulties stem from the possible “stiffness” of the reaction terms. Although many 
excellent numerical methods are now available for the nonreacting case (ip = 0) which give 
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high resolution and sharp shocks, it is not clear to what extent these methods can be used 
in the reacting case. 

The kinetics equations often include reactions with widely varying time scales. More- 
over, many of the chemical time scales may be orders of magnitude faster than the fluid 
dynamical time scales. This can lead to problems of stiffness akin to the classical stiffness 
problems of ordinary differential equations (ODEs). Stiff ODEs arise, for example, in mod- 
eling chemical kinetics in a uniform stirred reactor where the fluid dynamics terms drop 
out. The numerical difficulty with such problems is that some time scales will typically be 
much faster than the scale on which the solution is evolving and on which one would like 
to compute. This occurs when the fast reactions are in near-equilibrium during most of the 
computation. With many numerical methods, including all explicit methods, taking a time 
step appropriate for the slower scale of interest can result in violent numerical instability 
caused by the faster scales. 

Of course if the fast reactions are always in equilibrium it may be possible to eliminate 
these reactions from the system. In the extreme case one obtains equilibrium gas dynamics 
in which the kinetics are not explicitly modeled but the equation of state varies with the 
mixture. In many problems, however, nonequilibrium effects play an important role and 
must be included. 

In practice, it is common to take time steps which resolve the fastest scales. Boris and 
Oran[2] suggest as a rule of thumb that time steps must be restricted so that the energy 
release from chemical reactions does not change the total energy in any cell by more than 
10-20%. It is more desirable, however, to develop robust methods that can allow larger 
time steps. This will naturally incur some loss of resolution - reaction fronts will lose their 
structure and approach discontinuities, for example. What we should demand is that the 
correct discontinuities are obtained. They may be smeared out due to numerical diffusion, 
but should represent the correct jumps in the correct locations. This goal can be achieved 
for the nonreacting case by using “conservative” numerical methods[12] (See Section 4). 

In this paper we investigate the extent to which this goal can be achieved for equation 
(1) using various popular finite difference techniques. In particular, we introduce and study 
a simple one-dimensional scalar model equation which illuminates some of the difficulties 
sure to be encountered also in solving more realistic equations. We investigate the following 
questions: i) Can we develop stable methods? ii) Can we obtain “high resolution” results, 
with sharp discontinuities and second order accuracy in smooth regions, and iii) Do we 
obtain the correct jumps in the correct locations? 

Numerical stability is typically not a problem. A variety of excellent implicit methods 
have been developed for solving stiff systems of ODEs, and many of the same techniques 
can be applied to the stiff source terms in (1) to obtain stable methods for solving this 
system. 

The second question is investigated in Sections 2 and 3, where we will see that with 
some care, second order accuracy and reasonably sharp discontinuities can be obtained. 

The third question is the most interesting. For stiff* reactions it is possible to obtain 
stable solutions that look reasonable and yet are completely wrong, because the disconti- 
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nuities are in the wrong locations. Stiff reaction waves move at nonphysical wave speeds, 
often at the rate of one grid cell per time step regardless of their proper speed. 

This phenomenon has also been observed by Colella, Majda, and Roytburd[5] who made 
a similar study of the limiting behavior with increasing stiffness for various model systems. 
In particular they look at the Euler equations coupled with a single chemistry variable 
representing the mass fraction of unburnt gas in a detonation wave. These waves have the 
structure of a fluid dynamic shock that raises the pressure to some peak value, followed 
immediately by a reaction zone that brings the pressure back down to a new equilibrium 
value. On coarse grids it is not possible to resolve this combustion spike and the best one 
can hope for is a single discontinuity linking the two equilibrium values and moving at the 
correct speed. 

Colella, Majda, and Roytburd apply Godunov’s method and a high resolution extension 
of Godunov’s method[6] to this problem. The source terms are handled by splitting and 
solving the resulting ODEs exactly, so that stability is not a problem. However, they observe 
that on coarse grids the numerical solution is qualitatively incorrect. The computed solution 
consists of a weak detonation wave, in which all the chemical energy is released, followed 
by a fluid dynamic shock traveling more slowly. The reaction wave always travels at the 
speed of one mesh cell per time step, which is totally nonphysical. 

A simpler model system is also studied in [5] and is shown to exhibit similar behavior 
numerically. This system is essentially Burgers’ equation coupled with a single reaction 
equation. By studying this system and its numerical solution theoretically, progress has 
been made in understanding the structure of numerical solutions of the reacting Euler 
equations. 

However, the essential numerical difficulty can be identified and studied most easily by 
looking at even simpler equations. This same numerical behavior of discontinuities traveling 
at incorrect speeds can be observed in scalar problems. We have found it illuminating to 
study the model problem 

(2) u t + u x = tp(u) 


with 

(3) 4>(u) = -nu{u - l)(u - i). 

This is the linear advection equation with a source term that is stiff for large \i. Along the 
characteristic x = xq + t, the solution to (2) evolves according to the ODE 

(4) ^-u(x 0 + t,t) = ip(u(x 0 + t,t)) 

at 

with initial data u(a?o,0). This equation has stable equilibria at u = 0 and u = 1 and an 
unstable equilibrium at u = 1/2. For large [i and arbitrary initial data the ODE solution 
consists of a rapid transient with u approaching 0 (if u(x o, 0) < 1/2) or 1 (if u(x o, 0) > 1/2). 
Consequently, the solution u(x,t) to (2) with initial data w(x,0) rapidly approaches a 
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piecewise constant traveling wave solution w(x — t), where 


w(x) 


0 if u(z,0) < | 

< j if u(x, 0) = j 
, 1 if u(x, 0) > |. 


In particular, the solution with piecewise constant initial data 


( 5 ) 


«(x,0) 


{ 1 if x < xo 
0 if x > xo 


is simply u(x,t) = u(x - t, 0). In this case the ODE solution is in equilibrium on each side 
of the discontinuity, which theoretically behaves as it would if the source term were not 
present and we simply solved the linear advection equation u* + u x = 0. 

This linear discontinuity could easily be converted to a shock by replacing u x in (2) by 
f(u) x for some nonlinear flux function /. However, the numerical behavior is qualitatively 
the same in either case and nonlinearity of the flux is not the source of the difficulties of 
primary interest here. 

Other functions ip(u) could also be considered. A model corresponding more closely to 
the “ignition temperature kinetics” of [5] is obtained by using 

«„)=( ifu> | 

y ’ \ 0 if u < 


This gives numerical behavior similar to what is reported here for (3) and will not be 
considered further. 

Clearly this scalar model is inadequate as a full test of any numerical method. However, 
it does model one essential difficulty encountered in reacting flow problems, and is sufficient 
to point out difficulties that may arise also on more complicated systems of equations. 
Moreover, due to the simplicity of this equation, numerical problems that do arise can 
be easily understood and their source identified, yielding insight that may be valuable in 
developing better methods. 

We will discuss two different approaches to constructing numerical methods for (1) and 
compare their numerical behavior on the model problem (2) for various values of //. For 
simplicity we only discuss the one dimensional version of (1), in which < 7 = 0 , but in each 
case two-dimensional analogues are easily defined. Forward and backward differences of 
^-fluxes can be included in the predictor-corrector methods along with differences of the 
/-fluxes. 

The first method we consider is based on MacCormack’s predictor-corrector method 
for conservation laws[14]. This second order accurate method can be modified to include 
the source terms, which appear in each step of the method. Stiff source terms are usually 
handled in a semi-implicit manner to obtain stability with reasonable time steps. We 
have found, however, that one very natural and commonly used modification does not 
preserve the second order accuracy of the semi-implicit method on time-dependent problems, 
although steady states are accurately computed. Based on a truncation error analysis, we 
show how this can be easily rectified in Section 2. 
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In order to avoid oscillations near discontinuities, MacCormack’s method can be mod- 
ified by adding a flux-correction step motivated by the theory of TVD methods[20], [21]. 
We will compare two different forms of this correction. 

The second approach we study is the splitting method, in which one alternates between 
solving the conservation laws (with no source terms) in one step, and the stiff systems of 
ODEs modeling the chemistry (with no fluid motion) in the second step. This approach has 
certain advantages, in that high quality numerical methods exist for each of the subproblems. 
Combining these via splitting can yield stable, second order accurate methods for the full 
problem. This is demonstrated in Section 3. 

Numerical tests on the model problem (2) reveal that methods can be devised by 
either of these approaches that will be stable and second order accurate as the mesh is 
refined. However, for realistic choices of grid and time step, stiff reaction waves will have 
the nonphysical behavior described above. This is investigated in Section 4. 

2. Extensions of MacCormack’s method. MacCormack’s method for a system of 
conservation laws is a two step predictor-corrector method in which backward differences 
are used in the first step and forward differences in the second step (or vice versa). The 
method is easily modified to include source terms in an explicit manner and maintain second 
order accuracy[19]. For the one dimensional system 

(6) u t + f(u) x = ip(u) 


this explicit method takes the form 


( 7 ) 


At/f = ~\(W) - /((/;.,)) + WCJ) 
uf ] = Uf + AuV 

At/f = -|(/(^i) - f(uj 1} )) + ki>(uV) 
U? +1 = Uf + ^(At/| x) + Atfj 2) ). 


Here h is the grid spacing in x and k is the time step. Computing the truncation error for 
this method shows that it is second order accurate in both space and time, as the grid is 
refined with k/h fixed. 

Note that if we set f{u) = 0, so that (6) reduces to a system of ODEs, then (7) 
reduces to the standard two-stage Runge-Kutta method. Clearly this explicit method will 
be inadequate if the system is stiff, in that the time step k required for stability will be 
much smaller than desirable for accuracy. 

It is natural to try to improve the stability of the method by making it semi-implicit, so 
that the source terms are handled implicitly while the flux terms are still explicit. In order 
to avoid solving nonlinear systems of equations in each step, a linearly implicit method is 
frequently used. Methods of this form have been used by many workers (e.g. Bussing and 
Murman[3], Drummond, Rogers, and Hussaini[7], and Yee and Shinn[21]). This method 
takes the form 
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= U i 

Ui 


n 

e n 



ratio e n _!/e„ 

n 

5.449(-3) 

- 

3.617(-3) 


Kg™ 

1.514(-3) 

3.60 

6.383(-4) 


I 1 


3.88 



FH 

9.826(-5) 

3.98 

1.467(-4) 

1.71 


Table 1 

Computed errors with different choices of Uj 


( 8 ) 


I-\W{U?) 

I- 


A Uf 

uf) 

AC/f 


U? + 1 


1)) + w;) 

Uf + AcH 1} 

- wP)) + wfy 

U? + i(AC/j 1} + A£/j 2) ). 


The values of Uj and Uj are still unspecified. In most of the papers referenced above, 
Uj = Uj = uf ] is used, as motivated by (7). However, another possibility is to use 

Uj = Uj = UJ. A truncation error analysis for the method shows that this latter choice is 
in fact preferable, since it gives a method that is second order accurate in both space and 
time. The traditional choice is second order in space but only first order in time. Thus it 
gives second order steady state solutions but would only be first order in time for unsteady 
problems. 

Actually, in order to achieve second order accuracy overall, it is only necessary to use 
Uj = 17”. The choice of Uj is immaterial so long as i>'{Uj) — Uj l ) + O(k). In particular 

Uj = f/j 1 or Uj = Uj 1 ^ are both allowed. In view of this it seems most efficient to take 

Uj = U™ in practice, since then the matrix [I — \k^'(Uj)] need only be computed and 

factored once and the resulting factorization used in each fractional step. 

These statements are justified in the appendix, where we study the truncation error for 
this class of methods. To confirm this analysis numerically, we compute solutions to the 
model problem (2), (3) with n = 1 and initial data 

u(z,0) = i - tan -1 (10(a; - 0.3))/7t. 

We look at the computed results at time t = 0.3, when the solution is still smooth. Table 
1 gives the error e n obtained with each method (measured in the infinity norm) for various 
values of n (n = l/h and we use k/h = 0.75). We also list the ratios e n _i/e n , confirming 
our predictions of order of accuracy. 

Flux Limiters. The method (8) is spatially centered and hence will typically give 
oscillatory behavior on problems involving steep gradients. To minimize this problem, 
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one can introduce flux limiter terms into the method, as motivated by the theory of TVD 
methods. This is described in more detail in [20], [21] and so we will be brief in our description 
here. 

Let Aj +1 / 2 represent an average of f'(u) between ?7” and t/" +1 , e.g., the Roe approxi- 
mation [15]. Let Rj + 1/2 be the matrix of right eigenvectors of A J+1 / 2 and A J + 1 / 2 the vector 
of corresponding eigenvalues. Also, let 


, - V?). 

The components of this vector give the coefficients of the decomposition of the jump f/” +1 - 
Uj into eigenvectors of A J+1 / 2 - Corresponding to the /th component of this vector, a^ +1 / 2 > 
define a limited version by 


Q‘j+ 1/2 = minmod («j-l/2> a i+l/2 5 «j+3/2)- 


Other versions of this limiter can also be used (see [20], [21]) but we will restrict our attention 
to this one. Finally, we define 


/ _ k . j 

u 3+1/2 - l A j+1/2 


and the smoothed absolute value 

, x _ f \z\ if | 2 | > € 

q{Z) ~\ (z 2 + e 2 )/2e if \z\ < e 


for some positive parameter e. In our present example this plays no role, and in fact q(z) 
simply reduces to the absolute value in the context where we use it. 

We now replace the last line of (8) by 

Uf = Uf + |(A + A Uf ] ) 

and then set 


(9) c/; +1 = ^j 2) + [R J+ i/24> j+ i/2 - Rj- 1 / 2 ^- 1 / 2 } 

where <f>j+i /2 has components 

^'+ 1/2 = 1 / 2 ) “ (^j+i/ 2 ) 2 K a j+i /2 - Q‘j+ 1 / 2 )- 

Note that for smooth solutions, <*j +1 / 2 = 0(h) and «j +1 / 2 — Q l j+ 1 / 2 = 0(h 2 ) and is also 

smooth. The perturbation to U in (9) will then be 0(h 3 ), leaving the method second 
order accurate. Near discontinuities, however, this modification serves to introduce an 
upwind bias, dropping the method to first order accuracy but reducing oscillations. For 
a scalar problem with no source terms, the resulting method is TVD. When source terms 
are present, the true solution may no longer be TVD, and it is not clear what the correct 
theoretical criterion should be. 
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The above flux correction procedure can be modified by basing the correction terms on 
U W rather than U n . For example, in place of 1 / 2 we would use 

— ( i>( 2 ) — /rC 2 )^ 

This approach is advocated in [21] and has the advantage that U n need not be saved for 
this correction. However, this is no real advantage if we intend to save U n for the second 
stage of MacCormack’s method, as we have argued that one should do for time dependent 
problems. When based on rather than f/ n , the method is no longer strictly TVD on 
scalar problems without source terms. This might argue for the use of U n . When source 
terms are present it is not clear which approach is superior. Good results for a steady state 
reacting flow problem were achieved in [21] with corrections based on U ^ . This approach 
has also been successfully used for unsteady problems in the nonreacting case[l], [16], [23]. 

The experiments below indicate that for the model problem (2), limiting based on U n 
is preferable for small values of kfi but that limiting based on U ^ may be more robust for 
larger values of kfi. 

Numerical results on discontinuous data. If we carried the numerical computa- 
tions of Table 1 out to larger /, the solution would become nearly discontinuous. For larger 
values of fi this sharpening happens more quickly. To investigate the ability of this method 
to deal with propagating discontinuities we consider the following initial data for equation 
( 2 ): 


( 10 ) 


u(x,0) 


1 if x < 0.3 

0 if x > 0.3. 


We now use Uj = Uj = Uf exclusively and compare the effects of the different limiters (no 
limiter, limiting based on J7 n , limiting based on U^). We take h = 0.02, k/h = 0.75 and 
various values of fi. Note that due to scaling properties of the equation and method, results 
at time T with a particular value of fi can equally well be regarded as results with fi replaced 
by /z//?, for arbitrary /?, at time j3T with time step f3k and grid spacing f3h (with x rescaled 
so that [0,1] becomes [0,/?]). Indeed, the critical dimensionless parameters that determine 
the performance of the method are the mesh ratios A = k/h and the product kfi of the time 
step and reaction rate. The value kfi determines the stiffness of the system. When kfi is 
large, relaxation to equilibrium occurs on a time scale that cannot be temporally resolved 
on the grid. 

Figure 1 shows computed results at t = 0.3 for /z = 1,10,100, and 1000 (kfi = 
0.015,0.15,1.5, and 15). Each row of figures illustrates a different value of fc/z. The three 
figures in each row correspond to different choices of limiter. We see several interesting 
things from these graphs: 

• For small kfi (0.015) oscillations are visible if no limiter is used and to a lesser 
extent if the limiter is based on U^ 2 \ while limiting on U n gives monotone results. 
This agrees with what is expected for the pure convection case (kfi = 0). 

• For larger kfi (0.15 - 1.5), there is a slight overshoot in all cases, of similar mag- 
nitude regardless of the limiter. Note that for the case of no limiter there is less 


8 



QIO'O = <jro " 7/ ^ 
































oscillation here than with smaller due to the stabilizing effect of the source 
terms that tend to restore u towards 1. 

• For large kfi (15), limiting on U n appears to be unstable (there are large scale 
oscillations near x = 0.3 not visible in the figure) whereas limiting on U ^ or using 
no limiter gives stable results. In each case, however, the solution is completely 
wrong! The discontinuity has remained at its initial location x = 0.3 rather than 
propagating. 

Note that for kfi = 1.5 there is also some discrepancy in the location of the discontinuity. 
The speed of propagation is slightly too small. For intermediate values of kfi it is possible 
to obtain results with the discontinuity anywhere between 0.3 and 0.6. This phenomenon 
of wrong propagation speeds for large kfi will be discussed in more detail in Section 4. 

3. Splitting methods. The semi-implicit predictor-corrector method (8) attempts 
to handle the fluid dynamics and chemistry simultaneously. An alternative approach is to 
employ a time-splitting in which one alternates between solving a system of conservation 
laws, with no source terms, and a system of ordinary differential equations modeling the 
chemistry. In the simplest case this splitting takes the form 

(11) U n+1 = S^(k)Sj(k)U n . 

Here Sj (k) represents the numerical solution operator for the system of conservation laws 

u t + f(u) x = 0 

over a time step of length k, and S^(k) is the numerical solution operator for the ODE 
system 


u t = u ). 


To maintain second order accuracy, the Strang splitting[17] can be used, in which the 
solution J7 n+1 is computed from U n by 

(12) t/ n+1 = Sip(k/2)Sf(k)Sip(k/2)U n . 

Naturally, when several time steps are taken the adjacent operators 5^(A:/2) can be com- 
bined to give 

u n = s^kMSfikKs^Sfivr-'s^kwu 0 . 

In this form the method is nearly as efficient as (11). 

The splitting approach has also frequently been used to solve reacting flow problems[2], 
[4], [5]. At first glance it may appear to be less satisfactory than an unsplit method such 
as (8), since in reality the fluid dynamics and chemistry are strongly coupled and cannot 
be separated. However, the fact that the splitting (12) is second order accurate suggests 
that the interaction of different effects is adequately modeled by a split method, at least 
for smooth solutions. Moreover, there are distinct advantages to the splitting from the 
standpoint of algorithm design. High quality numerical methods have been developed both 
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for systems of conservation laws and for stiff systems of ordinary differential equations. By 
decomposing the problem into subproblems of these types, it is possible to take advantage 
of these methods directly. To some extent the mathematical theory that supports them 
can also be carried over. By alternating between using a high resolution method for the 
conservation law and a stable stiff solver for the system of ODEs, one can easily derive a 
method with excellent prospects of stability on the full problem. By contrast, attempting to 
devise a good hybrid method handling both effects simultaneously with good accuracy and 
stability properties can be difficult, as has been seen in the previous section. (But in the 
stiff case, we will still see the problem of incorrect wave speeds with the splitting method.) 

A split version of the method studied in Section 2 might take the form 


$*(*/ 2) : 


[l--kW?) J 

u* = u j 1 + a u; 


a v* 




S/(k) : 


(13) 


S^k/2) : 


w} 11 = - W-i)) 

uj 1} = u; + A{/j 1) 

Atf 1 = - /(tf >)) 

t/' 2) = U] + i(A + Atrf) 

uj- = uj 2 > + (*- +V! W - 


I- -ki>'(U")^ 

V n+1 _ V m* + AU * 


Am* = -k^m*) 


Here 4>* involves limited fluxes as before, based on U*. Alternatively, we can compute the 
limited value UJ* based on and replace R* and <f>* by and <f>^ respectively. 

Each of these methods could be replaced by other well known methods for the respective 
problems. For example, any implicit stiff solver, such as the trapezoidal method, could be 
used for S ^ and any of a wide variety of high resolution methods used for Sj. We consider 
the present form first as the logical choice for comparison with the previous results. The 
ODE method used in (13) for will be referred to as the “linearized implicit method”. 

Figure 2 shows the same set of experiments as in Figure 1, now with the splitting 
method. We observe that 

• For small kfi either choice of limiter (based on or U *) works well, and good 
results are obtained. 

• For k^i — 15 the discontinuity again moves at the wrong speed, now too fast. In 
fact, it has moved to x = 0.7 and so is moving at speed 4/3 rather than 1. Since 
k/h = 3/4, this indicates that the wave is moving at the speed of one mesh cell per 
time step. 

• A large overshoot occurs in one mesh cell behind the discontinuity for kfj, = 15, 
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/■ 7 ‘ = 0.015 


no limiter 


limit based on U' 


limit based on 





Fig. 2. Numerical results using splitting method with discontinuous initial data . 
computed solution. 


true solution , + : 
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regardless of the limiter used. 

With regard to this last observation, it appears that the overshoot must originate 
within the ODE-solving step. The flux-limiter method is applied only to the homogeneous 
conservation law and should give no overshoots, at least in the case where we limit based 
on U *. In other words, Sj(k) keeps monotone data monotone and therefore the lack of 
monotoniciy must be generated by S^(k). Note that this solution operator works pointwise 
(for example U* is a function only of t/", independent of U? for i ^ j), and so is oblivious 
to the gradient in u. What it does see, however, is a nonequilibrium value of u near 
the discontinuity. The linearized ODE method used in (13) is stable but converges in an 
oscillatory manner to the steady state of a stiff equation and we are seeing this here. In 
ODE terminology, the Sj, is not an L-stable method (see e.g. [8], [11]). 

These overshoots can be avoided by switching to a different ODE method. For example, 
if we leave Sj(k) unchanged but change S^(k/ 2) to the trapezoidal method, then these 
overshoots disappear for this value of kfi (Figure 3), but note that the propagation speed is 
still wrong. With the trapezoidal method we compute, for example, U* from Uf by solving 
the nonlinear equation 

Vj=U?+ + 

Although we obtain monotone profiles in Figure 3, the trapezoidal method also experiences 
overshoots if we go to still larger values of kfi. The use of an L-stable method such as the 
backward Euler method might eliminate this problem more generally, but backward Euler 
is only first order accurate. One might consider the use of higher order BDF methods (the 
“stiffly stable methods” of [8]), but the second order BDF method is already a two step 
method and in the present context we appear to require a one-step method, because of 
the nature of the splitting method. Implicit Runge-Kutta methods are a possibility. For 
reaction equations a special asymptotic method has been developed by Young and Boris[22] 
(see also [2]) which may avoid this problem. Another possibility is to use several steps 
of an ODE solver for S^(k/ 2), i.e., subdivide the time interval to a point where we can 
more adequately resolve the transient approach to equilibrium. It appears that this fails to 
achieve our goal of using time step large relative to the fast time scales, but note that we 
would need to do this refinement in time only in regions where nonequlibrium conditions 
hold. At grid points where u starts out close to equilibrium (e.g. those for which Ar^(u) is 
small, presumably most grid points), a single step of the linearized implicit method used in 
(13) is adequate to maintain stability. 

This is another advantage of the splitting method - since the ODE solver is decoupled 
from the fluid solver and is applied at each grid point independently, it is easy to change the 
ODE solver or even to use different solvers at different points depending on the character of 
the flow. This approach is also advocated by Young and Boris[22], who suggest using their 
asymptotic integration method at stiff points and explicit Euler elsewhere. 

We stress again, however, that improvements to the ODE solver cannot cure the problem 
of incorrect propagation speeds. In the next section we will investigate the source of this 
difficulty. 
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Fig. 3. Results for kp = 15 when the trapezoidal method is used. — ; true solution , ; computed solution. 


4. Nonphysical wave speeds. The numerical results presented above indicate a dis- 
turbing feature of this problem - it is possible to obtain perfectly reasonable results that 
are stable and free of oscillations and yet are completely incorrect. Needless to say, this 
can be misleading. In order to understand how the phenomenon occurs, it is sufficient to 
consider a simpler version of the splitting method, in which we use the splitting (11) with 
the first order upwind difference method for S / and the exact solution operator of the 
ODE for S^,. The method is then 

(H) u; = u? - t(i/" - i/jLj) 

Uj +1 = S^k)U*. 

We use the exact solution operator for Sj, to avoid the suspicion that difficulties are caused 
by the ODE solver. 

We also want to stress that for this scalar problem the splitting itself should not be 
suspect. In fact, it can be argued that the splitting method is the correct approach in the 
following sense. If 5/ and represent the exact solution operators, then they commute 
and the true solution u(x,t ) in fact satisfies 

u(x,t) = S^,(t)Sf(t)u(x, 0). 

This is just a restatement of the fact that the solution is obtained by integrating the ODE 
along characteristics, since Sj(t)u(x, 0) = u(x - f,0). But this says that the true solution 
at time t + k can be obtained from the solution at time t via the split method 

u*(x) = Sf(k)u(x,t) [=«(*-&,<)] 
u(x,t-\-k ) = S^,(k)u*(x). 


The method (14) is a direct discretization of this in which we replace 5/ by the upwind 
method. This amounts to replacing u(x - k, t) by the interpolated value 

(‘-sH+e 0 *- 

To see why this apparently reasonable method gives incorrect results when kfi is large, 
suppose we take initial data 


( 15 ) 



1 if j < J 
0 if j > J 
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for some J. Applying the first step of (14) gives 


(16) 



if j < J 
if j = J 
if j > J . 


5 


where A = k/h. In the second step we solve the ODE, which gives £/" +1 = UJ for j / J. 
For j = J the value we obtain depends on A and the size of ky. The interesting case is 
when ky >• 1, so that U j is restored to near equilibrium at the end of the time step. The 
equilibrium value reached depends on A, which by (16) is the initial condition for the ODE 
at grid point J . If A < 1/2 then the solution rapidly decays to zero and so Uj +1 « 0. If 
A > 1/2 then the solution rapidly approaches 1, so Uj +1 « 1. We thus obtain the following 
results, depending on the mesh ratio A: 


If A < 1/2: 
If A > 1/2: 



if j < J 
i i j> J 

if j < J + 1 
if j > J + 1 


The same behavior occurs in each time step and so we obtain a wave moving with speed 0 
if A < 1/2 or with speed 1/A (i.e. one mesh cell per time step) if A > 1/2. (If A is very close 
to 1/2 this argument is not valid. In fact, it is easy to verify that if A = 1/2 the method 
gives the correct speed 1. However, this is an unlikely special case.) 

In general we obtain propagation at a nonphysical speed that is purely an artifact of the 
numerical method. The problem lies with the smearing of the discontinuity caused by the 
advection, which introduces a nonequilibrium state ( Uj = A) into the calculation. Unfortu- 
nately, any conservative shock capturing method for the conservation laws will necessarily 
introduce some smearing since the true shock location almost never coincides with a cell 
boundary. At least one point in a shock is necessary in order to represent a discontinuity 
within a cell. As soon as a nonequilibrium value is introduced in this manner, the source 
terms turn on and immediately restore equilibrium, thus shifting the discontinuity to a cell 
boundary. 

It is difficult to see how this problem can be avoided using standard finite difference 
methods of the type used here, short of increasing the resolution considerably so that ky 
is small. To see how small ky must be to obtain reasonable results, it is interesting to plot 
the observed wave speed as a function of ky for fixed k/h. Because of the form of the true 
solution it is natural to define the wave speed in time step n by 


wave speed = | ^ U ? 1 j ■ 

For large ky this is essentially equal to 0 (if A < 1/2) or 1/A = h/k (if A > 1/2) in each 
step. For smaller ky this varies with n in a regular but generally oscillatory manner. To 
compare wave speeds for various ky, we define an average wave speed by averaging this 
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Fig. 4. Average wave speed as function of kp for A = k/h fixed 


function over a fixed time interval to to J n , or equivalently as 


(17) 


average speed = 


(!„ - fo) 


E v 7 - E V “ 


Figure 4 shows this average speed as a function of kfi for several values of A. 

This can be equivalently viewed as giving results on a fixed grid as /i varies or as giving 
results with fixed /i as k and h are varied (with A fixed), i.e., as the grid is refined. The latter 
viewpoint is more relevant to the discussion here, and shows that a substantial refinement 
(e.g. k/j. < 1) is necessary to obtain reasonable results. 

In these calculations k/h was held fixed. One might also consider fixing h but taking k 
much smaller in an attempt to resolve the nonequilibrium effects. Figure 5 shows that this 
is unsuccessful. Here h is held fixed at 0.01 and for various values of n the speed is plotted 
as a function of A = k/h, as k is varied. As expected, the correct speed is obtained only at 
A = 1/2 and A = 1. Note in particular that letting k — ► 0 with h fixed is detrimental and 
that if hfj, > 1 there is little hope of obtaining the correct speed. These results indicate that 
spatial resolution is as important as temporal resolution. This is not surprising, since it is 
the smearing of the discontinuity that is the source of the difficulty, and the extent of the 
smearing depends on the spatial resolution. 

If we wish to solve such problems without refining the grid to the extent indicated 
above, we must consider alternatives to the uniform finite difference methods considered so 
far. We must find methods that are capable of essetially increasing the spatial resolution 
without excessive refinement of the overall grid. (One could of course use local refinement 
only near the reaction fronts, but a large degree of refinement may be required. Unless 
such refinement is required for other reasons, this also seems wasteful.) Front tracking is 
one possibility, in which the reaction fronts are replaced by sharp discontinuities that are 
explicitly tracked as the solution evolves. This would probably give the best results, but 
is quite complicated in multi-dimensional problems. It would be nice to develop methods 
that can deal with stiff reaction fronts more robustly without requiring explicit tracking. 
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Fig. 5. Average wave speed as function of X for h fixed 


It is illuminating to compare the present situation with that of a homogeneous system 
of conservation laws with no source terms. In the latter case, the use of a “conservative” 
numerical method (as defined below) guarantees that an isolated numerical shock of the 
type considered here must propagate at the correct speed. It may be smeared out over 
several mesh points, but the speed, defined in manner analogous to (17), must be correct 
by conservation. To see this, note that by defining the cell average 

_ 1 /‘ I J+l/2 . , , 

Uj = T / u{x,t n )dx 

and integrating the conservation law u t + f(u) x = 0 over [*_,•_ t/2i x j+i/2] x [*n,fn+i], we 
obtain 

, , If ftn+l [t n + 1 

(18) ti" + = u] - - |jf f(u{x j+ 1/2, i )) dt - f i /2 , t)) dt . 

Summing this expression over j gives cancellation of the flux terms so that we are left only 
with fluxes at the boundaries of our region. 

Our numerical values 17” are approximations to u”. A finite difference method is said 
to be conservative if it can be written in the conservation form 

(i») v?" ' = v? - |(J7 +1 / 3 - c;_ 1/2 ] 

where F^ ±1 j 2 are the numerical fluxes based on U at neighboring points, and kF™ ±1 ^ 2 
approximates the corresponding integral in (18). Summing (19) over j gives the same 
cancellation of fluxes as in the true solution. Provided the fluxes at the boundaries of 
the region are correct, we maintain the correct total sum in each time step and hence the 
correct speed in the simple case of piecewise constant data with a single discontinuity. Lax 
and WendrofF[12] have shown more generally that a convergent conservative method must 
converge to a weak solution of the conservation laws, and thus must give discontinuities in 
the correct locations. 
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If we now include a source term and integrate 


+ f{u) x = ip{u) 

as before, we obtain 

1 r ftn+l rin+l 

(20) u] +1 = f(u(x j+1/2 ,t))dt- f(u(xj_i/ 2 ,t))dt 

1 rtn + i rx 3 + 1/2 

+ — / / ip(u(x,t))dx dt. 

hJtn 

The new term appearing here does not undergo cancellation when we sum over j, and 
consequently it is important that this term is modeled accurately if we are to obtain the 
correct behavior. 

Note that for the model problem we are considering, where the true u is everywhere 
in equilibrium except at the discontinuity, we have ip(u) = 0 almost everywhere and so the 
integral of ip in (20) is zero. In the numerical methods we have been considering, this integral 
is approximated by something analogous to kip{ u"). This is a reasonable approximation if 
u is smooth, but very poor in the present context. We are replacing the average value of 
ip(u) (which should be zero) by ip evaluated at the average value of u (which may be far 
from zero). 

One possible approach toward deriving better numerical methods is to attempt to model 
the integral of ip in (20) more accurately than simply using kip(U™). One possiblity is 
to compute some approximation v(x) to u(x,t n ) based on the grid values Uf, and then 
integrate ip(v(x)). One way to obtain a local reconstruction is to use the “subcell resolution” 
approach of Harten[10]. This method was originally proposed as a way to obtain sharp 
contact discontinuities in nonreacting flows, but appears promising in our context as well. 
The idea is to construct a piecewise polynomial function based on the data t/j 1 that may 
have discontinuities within the cells. Smoothness criteria and conservation are used to 
locate the discontinuities. Harten has tested a version of his method on the model problem 
considered here and reports excellent results[9], as would be expected on this scalar problem 
with piecewise constant solution. It is not yet clear to what extent this approach can be 
extended to systems of equations, and ultimately to multidimensional problems. 

5. Conclusions. We have proposed a very simple scalar equation as a model problem 
for understanding the behavior of numerical methods on reacting flow problems. We have 
considered two classes of numerical methods for this problem: MacCormack style predictor- 
corrector methods and splitting methods. In either case it is possible to derive second order 
accurate methods that are stable even for very stiff problems. However, all of these methods 
are subject to another numerical difficulty in the stiff case — incorrect propagation speeds 
of discontinuities. We have shown that this results from a lack of spatial resolution in 
evaluating the source terms. A nonequilibrium value in the numerical representation of the 
discontinuity, when viewed as the average value of u over a large mesh cell, will cause the 
source terms to be activated over this entire cell in a nonphysical manner. In order to avoid 
this difficulty, it will be necessary to increase the resolution of the discontinuity, at least for 
the purpose of evaluating ip(u). One possibility is to use some form of mesh refinement or 
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shock tracking. A more appealing alternative is to attempt to model the integral of in 

(20) more accurately using, for example, subcell resolution. It is not yet clear to what extent 
these approaches are practical for multidimensional systems of equations. The development 
of new methods along these lines is the subject of ongoing research. 

6, Appendix: Accuracy of the method (8). The truncation error T(x,f) of (8) is 
defined by 

(21) T(x,t) = jr(u(x,t + k) - u(x,t)) - ^(Au^(x,t) + 
where 

Au (1) (x,t) = I- ^ki/>'(u(x,t )) j ^~^[f(u(x,t)) - f(u(x - h,t))] + kip(u(x,t))^ 

u^(x,t) = u(x,t) + Au^\x,t) 

and 

Here the choice of u and u correspond to the choice of U and U in (8), i.e. either u(x , J) or 
t ). To compute the order of accuracy we must expand in Taylor series and simplify the 
expression (21) for T(x,t). This is easiest to do if we first consider the choice U = U = U n . 
Then (21) becomes 

T(x,t ) = i(«(x,/ + k)-u(x,t))+ -^-[1- ^V’ , («0M))] -1 

(22) x {f(u(x, t )) - f(u(x - h, t )) + f(u^(x + h, t)) - f(u^(x, t )) 

— 2h'ifi(u(x,t))}. 

Using the approximation 

|/ - ifcV>'(«(x, *)) = I + i kip'(u(x , *)) + 0(k 2 ), 

we obtain 

u^\x,t) - u(x,t)+ I+^kip'(u)-\ ] +krl>(u ) ) 

= u + k[ip(u) - f(u ) x ] + 0(k 2 ) 

= u + ku t + 0{k 2 ) 

where u = u(#,tf). Consequently, 

f(u^\x,t)) = f(u ) + kf\u)u t + 0(k 2 ) 
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Moreover, the 0{k 2 ) terms here are smooth functions of x and so will cancel to 0(k 3 ) when 
we compute f(u^(x + h, t)) — f(u^\x,t)), giving 

f(u^\x + h , *)) - f(u^\x,t)) = {[/(«) + hf(u)x + -h 2 f(u) xx + • • •] 

■\-k[f'(u) + hf'{u) x u x + • • • ][ut + hu tx + •••] + 0(A: 2 )| 
-{f(u) + kf'(u)u t + 0(k 2 )} 

= hf{u) x + ^h 2 f(u) xx + hk(f(u) x u t + f\u)u tx ) + 0(k 3 ) 
= hf(u) x + ~ h 2 f(u) xx + hkf(u) xt + 0(k 3 ). 

We also have that 

f(u(x,t)) - f(u(x - h,t )) = hf(u) x - i h 2 f(u) xx + 0(k 3 ) 
and so (22) becomes 

TOM) = [u t +^ku u + 0(k 2 )^ + ^\l+\k^u) + 0(k 2 )j 
x { 2 hf(u) x + khf(u) tx - 2 hip(u) + 0(fc 2 )} 

= u t + \ku tt + (/(«)* - u )) + \k[i>Xu)f{u) x + f(u) tx - ip'(u)4>(u)] + 0{k 2 ) 

= 0{k 2 ) 

since u t = ip(u) - f(u) x and u t t = ip'(u)u t - f(u) tx = rp'{u)^{u) - ip'( u )f( u )x ~ f( u )tx- This 
shows that the method (8) is second order accurate provided we use U j = U j = U ™ . 

Now consider what happens if we use Uj = instead of Uj = Uj 1 . Then the term 
2 hi>(u(x, t )) in (22) will become 

h[V>(«(M)) + ^(tt^(M))] = 2 hip(u(x,t)) + hi/)'(u(x,t))Au^Xx,t) H 

= 2 htp(u(x, t)) + hkij}'(u(x, t))u t + 0(k 3 ). 

Since this factor is multiplied by 1/h in computing T(x, f), this will cause an O(k) change in 
the truncation error and hence a reduction to first order accuracy in general. But note that 
in computing a steady state, where u t = 0, this perturbation drops out and so Uj = U 
can be used in that case. 

To justify our other claim, that alternative values of Uj are allowed provided that 
il>'(Uj) = i>'{U^) + O(k), consider the effect that using a different Uj would have on T( x, t ). 
For analytical purposes, we can rewrite (8) in this case as 

= -£[/(P") - /(t'j'-i)] + W?*) 

uf = U? + At/] 1) 

= -f[ /(^ , ,)-/(tfi 1 ’))+W“) 

A Uj = [/-i^(P j )]’ 1 [/-|t*'(P?)] Atrf 
UJ* 1 = UJ + i(A£/j 1) + AUj). 
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The first three lines are identical to the method already analyzed, i.e., (8) with Uj = Uj = 
Uj 1 . But now we compute a modified increment A Uj and use this to update Uj rather than 
A Uj 2 \ Clearly the method remains second order accurate provided AUj = A f/j 2 ^ + 0(k 3 ). 
But this follows easily by Taylor series expansion of the definition of AUj since ip'(Uj) = 
V(Uf) + 0(k) and AUf = 0{k). 
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